Characterising Complexity in Solar Magnetogram Data using a 
Wavelet-based Segmentation Method 
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ABSTRACT 

The multifractal nature of solar photospheric magnetic structures are studied using the 2D wavelet 
transform modulus maxima (WTMM) method. This relies on computing partition functions from the 
wavelet transform skeleton defined by the WTMM method. This skeleton provides an adaptive space- 
scale partition of the fractal distribution under study, from which one can extract the multifractal singular- 
ity spectrum. We describe the implementation of a multiscale image processing segmentation procedure 
based on the partitioning of the WT skeleton which allows the disentangling of the information concerning 
the multifractal properties of active regions from the surrounding quiet-Sun field. The quiet Sun exhibits 
a average Holder exponent 0.75, with observed multifractal properties due to the supergranular struc- 
ture. On the other hand, active region multifractal spectra exhibit an average Holder exponent -0.38 
similar to those found when studying experimental data from turbulent flows. 

Subject headings: Sun: flares, Methods: statistical, data analysis, Techniques: image processing, Magnetic fields, 
Turbulence 



1. Introduction 

Since the late 70 's and the propagati on of fractal 
ideas throughout the scientific community dMandelbrotl 
1982b , there have been numerous applications of the 
concepts of scale invariance, self- similarity, long- 
range dependence in many areas of physics, chem- 
istry, biology, geology, me teorology, economy, so - 
cial and material sciences (Aharonv& Feder 198 



| 1990l iLea-Cox & Wandll993l: Iwilkinson et al.lll995 
Taqqu et al.1 Il995l) . Unfortunately D F and H are 



9; 
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Westll 1 9901: Ivicsek et al.l[l994l: IBunde & Havlinll 1 994J ; 
Wilkinson et al 1995; F amily et al.ll99 5: Frisc hl 19951: 



global quantities that do not account for the pos- 
sibility of point-to-point fluctuations of the scaling 
properties of a fractal object. The multifractal for- 
malism was introduced in the mid-eighties to provide 
a statistical description of the fluctuations of regu- 
larity of singular me asures that a r e found in chaotic 
dynamical systems dHalsey et al.1 Il986t ICollet et aL 



Arneodo et alTll995ak IBunde et al.1 l2002h~ Various 
methods were developed to quantify scale-invariance 
properties through the computation of the fractal di- 
mension Dp for self- similar object s or the roughness 
exponent H for s e lf-affi n e fracta l s (IMandelbrotll 1982c 
Peitgen & Saupel 1 19871: iFeder! Il988l: lArgouletal. 



Il987l: iRandl Il989l) or in modelling of the energy 
casca d ing process in turbulent flows ([M andelbrot 



1974 IPaladin & Vulpianil 1 19871: iMandelbrotl 1 19891: 
Meneveau & Sreenivasanl Il99ll) . Box-counting and 



correlation algorithms were successfully adapted to 
resolve multifractal scaling for isotropic self- similar 
fractals by comp utation of the generaliz e d fract al 
dimensions D q (iGrassberger & Procaccial Il983allbl: 



Grassberger et al.1 1 19881). As to self-affi ne fractals, 
Parisi and Frisch (IParisi & Frischl Il985l) proposed, 
for the analysis of fully-developed turbulence ve- 
locity data, an alternative multifractal description 
based on the investigation of the sca ling behavior 
of the so-called struct ure functions (IFrischl 1 19951: 
Monin&Yagio^r]ll975h : S p (l) =< (6vi) p >~ fa (p 
integer > 0), where Svi(x) = v(x + /) - v(x) ~ l h ^ is 
an increment of the recorded signal over a distance /. 
Then, after rei nterpreting the roughness exponent as a 
local quantity dParisi & Frischl 1 9 85I:IMuzv et al.lll991 , 



19941: lArneodo et al.lll995cl) H 5vi(x) ~ l h{x) (power-law 
behavior), the D(h) singularity spectrum is defined as 
the Hausdorff dimension of the set of points x where 
the local roughness (or Holder) exponent h(x) of v 
is h. In principle, D(h) can be attained by Legendre 
transformi ng the structure f unction scal i ng exponents 
f„ dParisi & Frischl 119851: iMuzv et al.1 Il99ll Il994 : 
lArneodo et al.l Il995cl). Unfor tunately, as noticed by 
dMuzv et alJll99llll993lll994 . both the box-counting 
and structure function methodology have intrinsic lim- 
itations and fail to fully characterize the corresponding 
singularity spectrum since only the strongest singular- 
ities are a priori amenable to these techniques. As 
such, both meth ods are limited in their application 
to real data sets (JMuzy et al.lll994l: iGeorgoulisI 12005 : 
IConlon et alfeOOJT" 



In previous work, Arneodo and collabo rators (|Muzy et al 

199llll993l 1 1994 lArneodo et all 1 1995cl) have shown 



that there exists a natural way of performing a uni- 
fied multifractal analysis of both singular measures 
and multi-affine functions, which consists in using the 
continuous wavelet transfor m ( Goupil laud et al.l 1984; 
Grossmann & Morletl[l984l: lMeverlll990l: baubechiesl 



19921: iMallatl Il998h . By using wavelets instead of 
boxes, one can take advantage of the freedom of the 
choice of these "generalized oscillating boxes" to get 
rid of possible smooth behavior that might either mask 
singularities or perturb the estimation of their strength 
h. The other fundamental advantage of using wavelets 
is that the skeleton defined by t he wavelet tran sform 
modulus maxima (WT MM) dMallat & Zhongl ll992c 
Mallat & Hwan gill 9921) . provides an adaptative space- 



scale partitioning from which one can extract the D(h) 
singularity spectrum via the scaling exponent r(q) of 
some partition functions defined fr om the WT skele- 
ton. The so-called WTMM method dMuzv et all 1991 . 



1993L 119941 : lArneodo et al.lll995cl) therefore gives ac- 



cess to the entire D(h) spectrum via the usual Legendre 
transform D(h) = min q (qh - r(q)). We refer the reader 



to iBacry et al.l d 1993k and I JarTardl (11997) for rigorous 
mathematical results and to iHentschell dl994|) for the 
theoretical treatment of random multifractal functions. 
Applications of the WTMM method to ID signals 
have alrea dy provided insight into a wide variety of 
problems d Arneodo et al.l 120021) . e.g. the validation 
of the log-normal cas cade phenomenology of full 



developed turbulen ce d Arneodo et al.l Il998al 1 1999 



b; 



Delour et al.l 120011) and of high -r esolution tempora l 
rainfall dVenugopal et~aD l2006allbl : iRoux et al.l l2009h . 



the characterization and the understa nding of long 



range c orrela ti ons in DNA sequ e nces (lArneodo et al. 



1995bl 1 19961: lAudit et al.1 l200ll 120021) . the demon- 



stration of the existence of a causal cascade of in- 
formation f rom large to small scales in financia l 
time series dArneodo et al.lll998bl: IMuzv et al.ll2000h . 
the use of the multifractal formalism to discrim- 
inate bejvveejihealthy and s ick heartbeat dynam- 
ics dlvanov et al.l Il99a Il999|) . the discovery of a 
Fibonacci structural order ing in ID c uts of diffu- 



sion limited ag gregates dArneodo et al.l 11992 



nt ed ag g 
Kuhn et al.ll 19941) and in hard X -rav emission from so- 
lar flares dMcAteer et al]l2007l) . The WTMM method 



has been generalized from ID to 2D with the spe- 
cific goal to achieve multifractal analysis of rough 
surfaces u with f ractal dimens i on D F anywhere be- 
tween 2 and 3 (|Arrault et al.l 1 19971: lArneodo et al.l 
20001: iDecoster et al.ll2000l) . The 2D WTMM method 



has been successfully applied to characterize the inter- 
mitten t nature of cloud structure from s atellit e im- 
ages dArneodo et al.l Il999al: IRoux et al.l l2QQQh and 
to assist in the diagnosis of breast tissue lesi ons in 
digitized mammograms (Kesten er et al.l 120011) . In 
astrophysics, this method was adapted and used to 
characterize the anisotropic structure of atomic hy- 
droge n gas (HI) in the Galatic disk dKhalil et al. 



2006b . From the analysis of very large mosaics taken 



from the Canadian Galatic Plane Survey dTaylor et al. 



2003|) . directional roughness exponents were intro- 
duced to show that the HI in the Galactic spiral 
arms has a scale-dependent anisotropic signature 
while the HI in the inter- spiral arm regions exhibits 
scale-independent anisotropy. Along that line, the 
2D WTMM method was further applied to char- 
acteri ze the space- scale n ature of anisotropi c struc- 
tures dSnowetal.ll2008allbl : iKhalil et al.ll2009h and to 
perform objective segmentation of image features of 



! The fractal dimension of a rough surface associated to the graph 
z = S(x,y), where S(x,y) represents the height of the surface at 



location (x,y), is a quantity Dp between 2 and 3. 



interest from noisy backgrounds (IKhalil et al.l 12007 : 



Caddie et alJbOQ^ iRoland et alJl2009b We refer the 
reader to lArneodo et al.l (|2003|) for a review of the 
2D WTMM methodology, from the theoretical con- 
cepts to experimental applications. Recently, the 
WTM M method has been furth er extended to 3D 
scal ar dKestener & A rneodo 2003) as well as 3D vec- 
tor dKestener & Arneodol l2004 12007k fields analysis 
and applied to 3D (velocity, vorticity, dissipation, en- 
strophy) numerical data issued from direct numerical 
simulations (DNS) of incompressible Navier-Stokes 
equations. Because it combines singular- value decom- 
position and multifractal description, the so-called ten- 
sorial wavele t transform modulus maxi ma method for 
vector fields dKestener & Arneodoll2004 l2007h looks 
very promising for future simultaneous multifractal 
and structural (vorticity sheets, vorticity filaments) 
analysis of turbulent flows. 

Our aim here is to exploit the ability of the 
WTMM method to study compound systems that dis- 
play some non-analyticity in their multifractal spec- 
tra as the signature of some phase transition be- 
tween two underlying scale invariant components with 
different multifractal properties (IBohr & Tell 119881: 



Muzv et al.lll994l: lArneodo et al]fl995d) . These two 



components can both have some physical significance 
as previously experienced when using the WTMM 
method to d etect vorticity fila ments in swirling tur- 
bulent flows dRouxet al.lll999b or microcalcifications 
from br east tissue backgroun d in digitized mammo - 
grams dKestener et al.l 1200 ll: lArneodo et al.l l2003k . 



One of these components can be noise that may cause 
drastic distortions in the returned multifractal spectra. 
In this work we will follow a wavelet-based strategy 
inspired from the one previously used in ID to de- 
tect replication origins and promoters as jumps (dis- 
continuiti es) in ID noisy skew profiles in mammalian 
enomes (Brodie of Bro die et al .1120051: iTouchon et al. 



2' 



20051: [Nicolay et al. 2007) and in 2D to perform an 
objective and automatic segmentation of chromo- 
some territories in fluo rescence micros c opy imaging 
of mouse cell nuclei dKhalil et al J 120071: ICaddle et al. 



2007b and of gold formation on vapodeposited thin 



gold films dRoland et alJl2009h . 

The purpose of the manuscript is to demonstrate 
the suitability and reliability of the WTMM method 
to propose a wavelet-based segmentation procedure 
adapted to solar magnetogram data. In section [2 the 
basics of the 2D WTMM method are presented. Its 
ability to disentangle the underlying scale invariant 



components of a compound system displaying a phase 
transition in its singularity spectra is discussed and a 
strategy of segmentation is implemented. Section [3] is 
devoted to a test application of the proposed segmen- 
tation procedure on a theoretical data set with known 
multifractal properties. In section [4] we report the re- 
sults obtained when using this wavelet-based segmen- 
tation method to separate active regions from quiet- 
Sun features in solar line-of- sight magnetogram data. 
Our conclusion and future directions are then given in 
Section 

2. Segmentation methodology of compound mul- 
tifractal systems using the 2D WTMM method 

2.1. Basics of the 2D WTMM method 

The main steps of the 2D WTMM method are pre- 
sented here. Details can be f ound in lArneodo et al. 
J2000b:lDecosteretalJ d2000h : lArneodo et all d2003b : 
Khaliletal.ld2006h . 



1. Computation of the 2D continuous wavelet 
transform of the input image function f(x) with 
analyzing wavelets defined as the partial deriva- 
tives of a smoothing Gaussian kernel 0: 



TV[/](b,a) = 



Tft If] = a~ 2 fd 2 x i/n(a-\x - b))/(x) ^ 
W] = ^ fd 2 x iff 2 (a-\x - b))/(x) > 



= V{7>[/](b,fl)} 
= V{0 b ,,*/}, 

where if/\ = d(f>/dx, fa = d(p/dy and 0(x) = 
exp(-|x| 2 /2). Eq. (1) amounts to define the 2D 
wavelet transform as the gradient vector of f(x) 
smoothed by a dilated version cp(a~ l x) of the 
Gaussian filter. 

2. For each scale a, extract the WTMM edges de- 
fined as the locations b where the WT modulus 
Al^[/](b, a) = |T^[/](b, a)\ is locally maximum 
in the direction of the WT vector T^[/](b, a). 
These WTMM points lie on connected maxima 
chains. Along each of these maxima chains, 
locate the local maxima called WTMMM for 
WTMM maxima. Note that the two ends of an 
open maxima chain are not allowed to be a pos- 
sible WTMMM location. 

3. Extract the WT skeleton which is the set of max- 
ima lines £ Xo obtained by connecting these WT- 
MMM from scale to scale, starting at location 



(1) 



x at smallest scale. Start at the smallest scale 
a>min ~ 1 pixels (minimum size of the support of 
the wavelet function) and link each WTMMM 
to their nearest neighbor found at the scale just 
above. Proceed iteratively from scale to scale 
up to the largest scale a max (limited by the im- 
age size and border effects in wavelet transform 
computations). It is important to recall here that 
these lines contain all the information about the 
local Holder regularity properties of the function 
/ under consideration and that along a maxima 
line Xx that points to xo in the limit a — > + , the 
wavelet transform modulus behaves as a power 
law: 



AV[/][A>(0)]~" 



M*o) 



(2) 



where h(xo) is the Holder exponent, i.e. the 
strength of the singularity of the function / at 
the point xo. 

4. From the WT skeleton compute the partition 
functions: 

Z(q,a)= J] [M*[f](xe£,a)] q , (3) 

£e£(a) 

which allows to define the r(q) scaling expo- 
nents as 



Z(q, a) ~ a 






+ . 



(4) 



One can optionnally compute the companion 
partition functions h(q,a) and D(q,a) and de- 
fine the corresponding scaling exponents when 

h(q, a) = ZxeX(a) In |AWK X > a)\ W+[f](q 9 £, a) - 
D(q, a) = Zxex(a) Wflfiiq, £, a) ln(^[/](^, X, a)) 
where 
WflfKq, X, a) = [AV[/](x, a)] 9 /Z(q, a) (7) 

5. Compute the r(q) spectrum by performing lin- 
ear regression fits of In Z(q, a) vs In a and finally 
compute the D(h) singularity spectrum by Leg- 
endre transforming r(q): 



Note that alternative approaches to the WTMM 
method have been developed using discrete w avelet 



bases including the recent use of wavelet leaders (IJaffard et al 
20061: IWendt & Abrvl l2007[ IWendt et all l2007h . We 
think that the continuous WT better suits our goal 
to provide a selective multifractal analysis of multi- 
component images via some objective segmentation 
of maxima lines in the WT skeleton. 

2.2. Adapting the 2D WTMM method to the seg- 
mentation of compound multifractal systems 

For simplicity, we will assume that the compound 
multifractal systems of interest here can be considered 
as the sum of two scale invariant components: 



/(x) = / 7 (x) + / 77 (x) 



(9) 



characterized by the singularity spectra D\h) and 
D n (h) respectively. Ideally we will further suppose 
that D\h) and D n (h) have non-overlapping support 
\h l h 1 1 n \h u h 11 1 = and that h 1 < h 11 

in min> n maxi ' ' W'min' "max* ^ aiiU Lilcu n max ^ n min 

meaning that / 7 (x) possesses stronger singularities 
than f n (x). In the limit a —> + , the partition func- 
tion Z(q, a) (Eq. ©) can be split into two parts: 

Z(q,a) = Z I (q,a)+Z II (q,a) = Ci(q)a rI(q) +Cn(q)a TlI(q \ 

(10) 
where Cj(q) and Cn(q) are pref actors that depend on 
q. Since h I max < h 1 ^ , it follows easily that in this limit, 
there exists a critical value q crit so that: 



r(q) 



r\q) 
r n (q) 



for 
for 



q > qcrit 
q < qcrit 



(11) 



^faere^pje, the r(q) spectrum has a non-analyticity at 
q C riu .when crossing this critical value, there is a tran- 

"sftion^n-om one scale invariant component to the other. 
As illustrated in Fig. [U when Legendre transforming 
Eq. (fTTl) . one gets the upper envelop of the Z) 7 (/z) and 
D n (h) spectra, a classic al result for entropy in equilib 



rium statistical phys ics (Bohr & Tel 1988; Muz y et al. 



1994l;lArneodo et all 1995a) . In that respect the classi 



D(h) - min [qh - r(q)] . 
q 



(8) 



Alternatively D(h) can be computed from the es- 
timate of the scaling exponents h(q) and D(q) in 
Eqs. (0) and © respectively. 



cal 2D WTMM method does not provide separate ac- 
cess to D\h) and D H (h). 

Our segmentation strategy consists in using the WT 
skeleton to discriminate the maxima lines £}(a) as- 
sociated with singularities of / 7 (x) and maxima lines 
£ n (a) associated with singularities of / 77 (x), over the 
range [a min ,a max ] of accessible scales. This can be 
done theoretically by comparing the power-law be- 
havior of At^[/](x e JL) along each maxima line of 



the WT skeleton (Eq. ([2])). From the local estimate 
of /z(x), one can expect to partition the WT skeleton 
into two sub- skeletons, one made of the maxima lines 
X}(a) and the other one made of maxima lines £ n (a). 
In practice, this partitioning will suffer from the finite 
range of scales available to the analysis and the de- 
sired segmentation will require special care as far as 
finite-size effects and statistical convergence issues are 
concerned. 

3. Application of the wavelet-based segmentation 
method to synthetic data 

We consider an academic example im age with two 
fractal components: the fractal Dragon (lDudall2007|) 
embedded into a noisy background generated from 
fractional Brownian noise with Holder exponent H = 
-0.7. The fractal Dragon is a self- similar fractal 
defined as the limit set of an iterated function sys- 
tem (the Lindenmayer system) of the same type as 
the one used to g enerate the Sierpin ski gasket or the 
Von Koch curve (IMandelbrotl Il982h . but the fractal 
Dragon has less obvious geometrical symmetries. Let 
us note that the fractal Dragon is space-filling, mean- 
ing its fractal dimension is 2, whereas its boundary 
has a fractal dimension known analytically D Dragon = 

T , 1+ V73-6 V87+ V73+6 V87 v , COQ ^ A T r 

log 2 ( — 1 — — — ) 2* 1.5236. A sample frac- 
tal Dragon is shown in Figure [2 a) whereas the corre- 
sponding noisy two-component image is shown in Fig- 
ure (3d). As previously discussed, the wavelet analy- 
sis proposed in this work is sensitive to singularities, 
i.e. to points in the images where the signal is sin- 
gular. We expect the WTMM analysis of the image 
shown in Fig. [2jd) to simultaneously reveal multifrac- 
tal information about both the boundary of the frac- 
tal Dragon and of the rough background texture. Let 
us recall that the two components have known mono- 
fractal type self- similar properties, i.e. a singularity 
spectrum degenerated to a single point: (h = -0.7, 
D - 2) for the fractional Brownian noise and (h = 0, 
D = D Dragon =* 1.5236) for the boundary of the frac- 
tal Dragon. The roughness H = -0.7 of the fractional 
Brownian noise was chosen to mimic the texture of the 
quiet-Sun images (see Section [47Tb . Figures Ob) and 
[2e) illustrate the results of the computation of the WT 
maxima chains at the smallest scale; the arrows corre- 
spond to the WT vectors (Eq. ®) at the WTMMM lo- 
cations. Figures [2c) and [20 show the maxima chains 
at a scale twice as large as in Figures Ob) and He). 
When going from large to small scales, whereas the 




(b) r D(h) 



Fig. 1. — Illustration of a phase transition in the 
multifractal spectra of a compound system (Eq. ©). 
(T 7 (g), D\h)) and {T n {q),D n {h)) are the multifractal 
spectra for / 7 (x) and f n (x) respectively. The 2D 
WTMM method and more generally the multifrac- 
tal formalism, give access to the dashed r(q) curve 
(Eq. (fTTt ) in (a) and via the Legendre transform 
(Eq. ([5])) to the dashed D(h) sprectrum in (b) which 
is the supremum of the D\h) and D u (h) spectra. 




Fig. 2.— (a) The fractal Dragon (1024 x 1024). (b) 
WTMM chains at the smallest scale a = crw - 7 pix- 
els in a small 100 x 100 region of the fractal Dragon, 
(c) Same as in (b) for scale a - 2cr w . (d) The frac- 
tal Dragon embedded in a fractionnal Brownian noise 
(H = -0.7) background of twice as large amplitude, 
(e) and (f) are the same as (b) and (c) but for the noisy 
fractal Dragon. In (b,c,e,f), the black arrows represent 
the WT vectors originating at the WTMMM. In (b) and 
(e) the background image is the smooth-convoluted 
image 0b,« * / at scale a = crw- 



boundary of the fractal Dragon is better and better ap- 
proximated by some WTMM chains (edge detection 
in the smoothed image), an increasing number of addi- 
tional maxima chains start emerging as the signature of 
the presence of a colored noise (Fig.Oe) as compared 
to Fig. Eft)). 

As previou sly emphasized dMuzy et al.l 19941; lArneodo etal 
1995ci 120031) . the set of maxima lines that defines 
the WT skeleton contains the space- scale information 
necessary to recover the underlying multifractal prop- 
erties. In Figures [2a) and[2b) are shown in a logarith- 
mic representation, the behavior of the WT modulus 
along the maxima lines computed for a noisy frac- 
tal Dragon with a noise amplitude respectively twice 
and five times as large as the fractal Dragon. Since 
the fractional Brownian noise is e verywhere singula r 
with Holder exponent H = -0.7 dMandelbrotlll982h . 
maxima lines pointing to noise features at small scale 
are characterized by a WTMMM power law behavior 
Mijf[f](a) ~ a~ 0J , while lines associated to the frac- 
tal Dragon boundary can be distinguished by the fact 
that Mijf[f](a) ~ a° ~ Const (no scale dependence). 
This leads us to implement the following segmenta- 
tion procedure: the space (log 2 a,log 2 M l /,[f](a)) is 
divided in two regions separated by a straight line of 
slope -0.7 < h s < and intercept log 2 M s . As shown 
in Figure [2a), for low noise amplitude, all the max- 
ima lines along which log 2 At^[/](a) decays slower 
than h s log 2 a when increasing a, are colored in red 
and associated with the Dragon boundary. On the con- 
trary, all the maxima lines along which log 2 M^lf^a) 
decays faster than h s log 2 a when increasing a, are 
colored in blue and associated with the noise com- 
ponent. But as shown in Figure [2b), for large noise 
amplitude the distinction of the two sub- skeletons is 
much more tricky at small scales where some entan- 
gling is observed. We thus adapt the segmentation 
criteria towards the largest scales in fully analogy 
with a different but conceptually simila r adaptation 
of the 2D WTMM segmentation method dKhalil et al. 



2007|) . Each maxima line is characterized by a length, 
i.e. its maximun scale a max and the WT modulus 
M$[f](a max ) at that scale. A given maxima line is said 
to belong to the Dragon sub-skeleton, if it satisfies the 
following condition : 

log 2 AV If] (a m ax) > h s log 2 a max + log 2 M s . (12) 

In Figures HJa) and H{b) are reported the results of 
the computation of the partition functions for the frac- 
tal Dragon alone and its noisy version after applying 





Fig. 3. — Log-log plot of WT modulus along the 
skeleton maxima lines versus scale. Lines are colored 
according to the segmentation procedure (Eq. (fT2l) ) : 
fractal Dragon boundary (red) and fractionnal Brown- 
ian noise (blue), (a) Noisy fractal Dragon with a noise 
amplitude twice as large as the fractal Dragon (see 
Fig. [2d)), (b) Noisy fractal Dragon with a noise am- 
plitude five times as large as the fractal Dragon. The 
dashed black line represents the segmentation condi- 
tion (Eq.[T2). 







-2024 




2 


: 


* 


4 

1 , - 


8 
6 


1.6 




4 


: 


;"■' , , , ° 



-0.8 -0.4 



Fig. 4. — Multifractal analysis of the fractal Dragon 
(o) and of the noisy (// = -0.7) fractal Dragon (•) 
after applying the segmentation procedure (Eq. (fT2l)). 
(a) h(q, a) vs log 2 a for different values of q\ the solid 
lines correspond to linear regression fits over the range 
of scales a e [2°,2 4 ] cr w (resp. [2 a5 ,2 35 ] cr w ) for 
the fractal Dragon (resp. the noisy fractal Dragon af- 
ter segmentation). The symbols (■) correspond to the 
h(q = 0, a) partition function obtained for the noisy 
fractal Dragon without any segmentation, (b) D(q, a) 
vs log 2 a. (c) r(q) vs q\ the dashed horizontal line is 
the theoretical spectrum r(q) - -1.5236 (Vg) of the 
boundary of the fractal Dragon, (d) D(h) vs h\ the sym- 
bols (•) correspond to the segmented D(h) spectrum 
for the fractal Dragon component; the (x) correspond 
to the extracted D(h) spectrum of the colored noisy 
background. The zoom-in window enlarges D(h) data 
corresponding to the noisy fractal Dragon (•) after ap- 
plying the segmentation procedure. 



the segmentation condition (Eq. (fT2l) ) with h s = -0.5 
and log 2 M s = -3.2. In Figure SJa), the partition func- 
tions h(q, a) (Eq. ©) of the noisy Dragon display a 
well defined scaling behavior over 3 octaves (com- 
pared to 4 octaves for the original fractal Dragon), for a 
wide range of values of q e [-2,3]. For negative q val- 
ues, at very small scales, the segmentation procedure 
fails to disentangle the two components due to the con- 
tribution from the noisy component with h(q) ^ -0.7. 
Nevertheless the gain in scaling is unquestionable as 
compared to the behavior of the h(q, a) partition func- 
tions without segmentation (the ■ in Fig.|4ta)). With- 
out any segmentation the partition functions are a mix- 
ture of different scaling behaviors from which reliable 
quantitative information cannot be extracted. In Fig- 
ures SJc) and SJd) are shown the corresponding r(q) 
and D(h) spectra. Despite some slight departure from 
monofractality for the segmented noisy fractal Dragon 
(that is also observed but to a lesser extent in the orig- 
inal fractal Dragon as the result of finite- size effects), 
one recovers a rather good estimate of the fractal di- 
mension D F = 1.57+0.03 of the fractal Dragon bound- 
ary. Furthermore, as reported in Figure EJd), our seg- 
mentation procedure has proved to be very efficient 
to estimate separately the D(h) singularity spectra of 
both the fractal Dragon and the noisy background. 
This efficiency is illustrated in Figure \5\ where a 3D 
(x, y, scale) space- scale visualization of the maxima 
chains of the noisy Dragon prior (Fig. [3a)) and after 
(Fig. Ob)) segmentation clearly confirms the elimina- 
tion of noise-induced small scale features that would 
otherwise severely affect the multifractal analysis. 





Fig. 5. — 3D visualization in the space-scale 
(x, v, scale) representation of the WTMM chains com- 
puted from the image shown in Fig. Od) before (a) 
and after (b) the segmentation procedure (Eq. (fT2l)). 
At each scale a, only the maxima chains containing 
at least one WTMMM belonging to the resulting WT 
skeleton are displayed. 



4. Application of the wavelet-based segmentation 
method to Solar magnetogram data 

Magnetic field measurements were obtained by the 
Michelson Doppler Imager (MDI) on the Solar and 
Heliospheric Observatory (SOHO), which images the 
Sun on a 1024x1024 pixel CCD ca mera through a 



series of increasingly narrow filters (IScherrer et al. 



19951) . The final elements, a pair of tunable Michel- 
son interferometers, enable MDI to record filtergrams 
with an FWHM bandwidth of 94 mA. In this paper, 
96-minute magnetograms of the full disc were used, 
which had a pixel size of ~2". For the purposes of 
this work, a series of magnetograms have been ana- 
lyzed to examine the difference in fractal properties 
between quiet and active solar regions. A total of 29 
magnetograms representative of the quiet Sun were 
taken from December 21 to December 22, 2006 and 
a similar series of 28 images representative of the ac- 
tive Sun were taken from October 27 to October 29, 
2003. In the solar photosphere, the large magnetic 
Reynolds number (~ 10 7 - 10 9 ) means that magnetic 
field lines will be ad vected with the flow of plasma 
(IMcAteer et alJ l2009|) . This system naturally leads 
to self- similarity, sugg esting a m ultifrac tal study is 
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appropriate ([Lawrence et al] 119931: lAbramenko et al. 
20021: lMcAteeretalJl2005h . As already mentioned, 



previous method of calculating the multifractal proper- 
ties of solar magnetic features are dependent on image 
resolution, thresholding, and instrument sensitivity. 
The WTMM method calculates the multifractal spec- 
trum of solar magnetic features based on the distribu- 
tion of gradients within the image at various scales. As 
such, the WTMM multifractal parameters are less sen- 
sitive to changes in image resolution and instruments 
than traditional methods. 

4.1. Quiet-Sun multifractal properties 

Examples of quiet and active MDI magnetograms 
analyzed are shown in Figure [6] (top left and top right 
respectively), with a histogram of the wavelet trans- 
form modulus Mip[f](b,a) at the smallest scale (bot- 
tom). Active regions result from an increased propor- 
tion of large magnetic elements of opposite polarity 
in close proximity to each other. The resulting neu- 
tral or magnetic inversion lines can be detected us 



ing s tandard wavelet-based techniques (llreland et al. 



20081) . As such active regions should contain a greater 
number of higher magnitude WT gradients. This is 
shown in Figure [6] (bottom), which suggests that mod- 



uli with values larger than 40 are unlikely in quiet- Sun 
magnetograms. Due to the different scaling properties 
of active regions and their surrounding quiet Sun, our 
goal is to segment the WT skeletons using the condi- 
tion defined in Eq. (fT2l) . As outlined in Section [2 and 
illustrated on synthetic data in Section [31 this should 
allow us to study the multifractal properties of active 
regions in a quantitative manner. 

The results of the computation of the multifractal 
spectra when averaging the partition functions over a 
set of 30 (505 x 505) quiet- Sun images without ap- 
plying the segmentation are reported in Figure [7] As 
shown in Figure |3a) and [2b), h(q, a) (Eq. ([5])) and 
D(q, a) (Eq. ©) display convincing scaling behavior 
over almost four octaves for q e [-2, 3] (symbol (o)). 
Linear regression fits of the data yield the non-linear 
r(q) spectrum shown in Figure Etc). This multifractal 
diagnosis can also be observed in Figure [T^a) where 
the slope h(q) of the partition function h(q, a) versus 
log 2 a definitely depends on q. The corresponding 
multifractal spectrum D(h) is shown in Figure [3d). 
From the top of the D(h) curve, we can see that quiet- 
Sun images are everywhere singular (D(q = 0) = 2) 
with a corresponding Holder exponent h(q = 0) ^ 
-0.75. The multifractality can be quantified by the 
so-called intermittency coefficient c^ that characterizes 
the width of the D(h) curve. As shown in Figures |3c) 
and Uld), the r(q) and D(h) data of the quiet-Sun im- 
ages are well fitted by a parabola 

r(q) = -c + Ciq- —q 2 , DQi) = c , 

2 2c 2 

(13) 
where Co - 2, c\ =* -0.75 and c^ - 0.22. Let us 
point out that quadratic multifractal spectra are pre- 
dicted by the so-called log-normal model that has been 
populari zed by the fu l ly-developed turbulen c e com- 
munity dFrischl 119951: lArneodo et all Il998al Il999bl: 
Delour et al.l l200l|) . In the present case, there is no 
particular evidence of the relevance of this model ex- 
cept that the observed r(q) and D(h) multifractal spec- 
tra are well characterized by their log-normal quadratic 
approximations. 

In order to understand the source of this intermit- 
tency, a upper threshold was imposed on each MDI 
magnetogram of the quiet Sun (Figure [8]). The thresh- 
old operation has the effect of removing large magnetic 
features resting on the boundary of the super-granular 
structures of the Sun. In Figure Etd), we can see that 
the multifractality of the thresholded quiet- Sun image 
(symbol ■) set is strongly reduced but not totally can- 
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Fig. 6. — (Top, Left) MDI magnetogram taken on 
December 20, 2006; (Top, Right) MDI magnetogram 
taken on October 28, 2003. (Bottom) Histogram val- 
ues of the wavelet transform modulus At^[/](a) at the 
smallest scale (cr w = 7 pixels) for MDI magnetogram 
images of a quiet Sun (solid) and active Sun (dashed). 
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Fig. 7. — Multifractal analysis of a set of 30 quiet-Sun 
images (505x505) prior (o) and after (■) thresholding 
(a sample thresholded magnetogram is shown in Fig- 
ure [8]). (a) log 2 h(q, a) vs log 2 a for different values 
of q\ the solid lines are linear regression fits over the 
range of scales a e [2°, 2 3J ] cr w . (b) log 2 D(q,a) vs 
log 2 a. (c) r(q) vs q\ the dashed straight line is the 
theoretical linear spectrum r(q) = -0.15q - 2 of the 
2D fractional Brownian noise with Holder exponent 
H = -0.75. (d) D(h) vs h; the dashed lines delimit the 
position of the top of the D(h) curves. Error bars cor- 
respond to standard deviation in the linear regression 
procedure. 



celled. We can also note that the average Holder ex- 
ponent is slightly shifted from < h >= c\ - -0.75 to 
-0.82, and the intermittency coefficient is reduced to 
C2 ~ 0.06 (Let us recall that a value c^ - means 
that the underlying process is monofractal). With- 
out the super-granular magnetic structure, the quiet- 
Sun multifractal spectrum looks much more monofrac- 
tal. This suggests that the magnetic features rest- 
ing on the boundaries of the super-granular structure 
are a major actor in the obs erved intermittent struc - 
tural properties of the Sun (iGeorgoulis et al.l 120021) . 
Since current models for the solar dynamo use in- 
format ion on the fractal d imension of solar disk as a 
whole (IPontieri et al.l2003|) . these new informations on 
the photosphere and the characteristic make-up of the 
quiet Sun should be incorporated in further theoretical 
works. 



4.2. Solar magnetogram active region segmenta- 
tion 

In this section, we highlight the use of the WTMM 
segmentation method on Solar magnetogram data with 
active regions, to demonstrate its ability to analyze the 
underlying multifractal properties of the active regions 
that are embedded in the surrounding quiet- Sun tex- 
ture. 

A sample 505 x 505 magnetogram MDI image con- 
taining an active region is shown in Figure |3 Fig- 
ures |3b) and [3c) show respectively the results of the 
computation of the WTMM chains before and after the 
segmentation at scale a - cr w ~ 7 pixels of a small 
150 x 150 excerpt focused on the active location. As 
explained in Section 14.11 WT skeleton maxima lines 




Fig. 8. — 256x256 quiet- Sun images. Image on the 
right is a threshholded version of the left one. Pixels 
with large absolute magnetic flux are shrinked down. 
Multifractal properties of quiet- Sun images and the 
corresponding thresholded versions are shown in Fig- 
ure [71 



associated with quiet-Sun structures have a character- 
istic scaling behavior described by M,j,[f](a) ~ a~ 0J . 
This behavior is used to derive the parameters char- 
acterizing the line h in Figure [TOl that will allow us 
to discriminate in the WT skeleton, the maxima lines 
(blue) that correspond to quiet- Sun structures: 

log 2 M^[f](a max ) < h Q log 2 a max + log 2 M Q , (14) 

where -0.7 < Hq < 0, so that for the selected maxima 
lines log 2 M^[f](a) will decrease fast enough when 
increasing the scale a to correspond to the quiet phase. 
To select the maxima lines (red) associated with the 
active region, we use another separating line l\ in Fig- 
ure na 

log 2 M#[f](amax) > h A log 2 a max + log 2 M A , (15) 

where this time < h A < 0.38, to limit the de- 
crease (if any) of log 2 M^lf] when increasing a. Note 
that the lines l\ and h have different slopes because 
some maxima lines cannot be clearly associated with 
either the quiet Sun or the active region state. In- 
deed those maxima lines are associated to features lo- 
cated near the boundary between quiet- Sun and ac- 
tive regions. When going from small scales to large 
scales the support of the analyzing wavelet starts cov- 
ering partly both regions preventing accurate classifi- 
cation. As illustrated in Figure [TTJ when fixing the 
segmentation parameters to Hq = -0.40, h A = 0.25 
and log 2 Mg = log 2 M,4 = 5.2, the space-scale na- 
ture of the methodology allows us to disentangle max- 
ima chains associated with the active region from those 
corresponding to the quiet Sun. Let us note that in a fu- 
ture work on large data sets, an automated parameters 
adjustment will be implemented using a clustering al- 
gorithm. In addition, a wrong choice of segmentation 
parameters can be observed in the partition function 




Fig. 9. — (a) 505x505 Active Region example image 
(October 28, 2003). (b) WTMM chains at the smallest 
scale (a = cr w pixels) in a small 150 x 150 region sur- 
rounding an active spot, (c) WTMM chains left after 
the segmentation step. 




Fig. 10. — Log-log plot of the WT modulus along 
the skeleton maxima lines versus scale. The dashed 
line denoted l\ with slope Ha = 0.25 and intercept 
log 2 Ma = 5.2 is used to identify WT skeleton max- 
ima lines associated to the active region (red). Ac- 
cording to Eq. ([T5K these lines have an ending point 
at highest scale (log 2 a max , log 2 M^[f](r, a max )) above 
l\. The dashed line denoted h with slope Hq = -0.40 
and intercept log 2 Mq = 5.2 is used to identify maxima 
lines associated to the quiet Sun (blue). According to 
Eq. (fT4l) . these lines have an ending point below h. All 
other lines (grey) are not clearly identified to belong 
either to the active site or the quiet surrounding. The 
values of the segmentation parameters Jia, log 2 Ma, Hq 
and log 2 Mq were chosen by examining the WTMM 
histogram at the smallest scale to extract at best a sub- 
skeleton specific to the active region. 



plots (not shown here) which display a phase transition 
phenomenon, i.e. a scaling behavior that changes from 
one state to the other when going from small scales to 
large scale due to to non-homogenous phases and mis- 
classified maxima lines. 

In Figure [12] are shown the results of the partition 
function computation for a set of 5 active region mag- 
netogram images taken on October 28 ?/l , 2003 (one 
image out of this set is shown in Figure Oa)). Parti- 
tion functions are computed separately for each sub- 
skeleton corresponding to the extracted action region 
maxima lines and quiet- Sun maxima lines shown in 
Figure \TU\ From these plots, one can see that the 
log 2 h(q, a) versus log 2 a data are nicely modeled with 
linear regression fits with slope h(q) that depends on q 
as the signature of multifractal scaling. This demon- 
strates that the segmentation procedure successfully 
extracts from the magnetogram images two scale in- 
variant components with different multifractal proper- 
ties. The r(q) (Fig.Ete)) and D(h) (Fig.[Ead)) spectra 
computed from the set of quiet-Sun maxima lines (blue 
lines in Fig. [10| are in good agreement with the ones 
previously computed in Section 14.11 from pure quiet- 
Sun images (Figs. [3c) and [3d) respectively). Us- 
ing the log-normal approximation (Eq. (fT3l)), we get 
Co = 2, c\ = -0.65 and c^ - 0.10. This means that the 
extracted quiet Sun appears (like the thresholded MDI 
magnetograms of quiet Sun in Figs |7fc) and Uld)) a 
little less intermittent as compared to the previous es- 
timate C2 = 0.22. As for the active region, the corre- 
sponding partition functions computed from the set of 
active phase maxima lines (red lines in Fig. [10| dis- 
play very convincing multifractal scaling behavior as 
quantified by the r(q) and D(h) spectra shown in Fig- 
ures Etc) and[l2td) respectively. Again these spectra 
are well approximated by the quadratic log-normal for- 
mula (Eq. (TT3l) ) with parameters Co = 2, c\ - 0.38 and 
C2 = 0.12. This indicates that the singularities associ- 
ated with the active region are space-filling (they are 
distributed on a set of fractal dimension D F = Co = 2), 
with a mean strength h(q = 0) = c\ = 0.38 meaning 
that magnetogram images can be considered as contin- 
uous on active regions (over the range of scales of our 
analysis) but noisy on quiet regions. 



5. Conclusions 

Many complex physical systems analyzed using 
fractal and multifractal techniques are surrounded or 
embedded in a noisy background sometimes originat- 
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Fig. 11. — 3D visualization in the space-scale 
(x, y, scale) representation of the WTMM chains com- 
puted from the image shown in Fig. [3a) before (left) 
and after (right) the segmentation procedure. The seg- 
mentation conditions are defined in Figure [lOl The 
maxima chains displayed contain at least one WT- 
MMM belonging to the resulting skeleton. 



ing from instrumental noise. As such these systems 
are a statistical combination of two distinct self- similar 
structures. This work addressed the need for an ac- 
curate calculation of the multifractal parameters of 
such complex systems. The presence of compound 
scale-invariant structures can result in an inaccurate or 
skewed calculation of the fractal and multifractal pa- 
rameters when studied as a whole. Using a wavelet- 
based multi- scale segmentation method, we show that 
it is possible to disentangle to some extent these two 
processes and accurately (up to finite-size effects) re- 
cover the multifractal characteristics of the system of 
interest. A theoretical test example for this method 
was provided in section [2 The removal of informa- 
tion relating to the background noise was highlighted 
in Figure [5] The quantitative results reported in Fig- 
ure 01 attest of the ability of this segmentation method 
to recover the multifractal parameters in question. 

Let us emphasize here that the application of this 
method to experimental data for which we do not have 
a priori knowledge of the possible underlying multi- 
fractal processes is a difficult task that requires much 
attention to perform the most objective segmentation 
which can not be infered or guided by some physi- 
cal rule or information. The multifractal analysis is 
a statistical tool that has direct connections with sig- 
nal and image processing, but not necessarily with the 
physics of the system per se. As noticed in section l4~2l 
the use of a clustering algorithm should greatly help 
in adjusting the multifractal parameters of the differ- 
ents components as well as in providing an automated 
procedure for processing large data sets. This will be 
reported in a future publication. 

The application of this wavelet-based methodology 
to quiet- Sun data has revealed the multifractal nature 
of this intermittent noisy component (< h >= -0.75) 
as mainly resulting from the super-granular magnetic 
structures. The quiet-Sun study was also necessary 
to get expertise for further analyzing more complex 
images that involve a segmentation before being able 
to clearly identify the underlying multifractal proper- 
ties. We have checked that the partition functions com- 
putations for the segmented quiet- Sun phase provide 
(i) convincing scaling properties and (ii) multifrac- 
tal spectra r(q) and D(h) estimates in good numerical 
agreement with the one measured in the previous cal- 
ibration step. The assumption of two non-overlapping 
D(h) is not inconsistent with the data. Let us notice 
that a phase transition can be observed in the partition 
function log-log plots when inappropriate segmenta- 
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tion parameters are chosen. In the case of overlapping 
D(h), there would be a large set of maxima lines in 
the WT skeleton that could not be genuinely sorted, 
which would prevent from building accurate multifrac- 
tal D(h) measurement. From the analyzed data, we 
were not able to distinguish more than two phases. 
Finally, this gives us good confidence in the segmen- 
tation proposed for solar magnetogram containing an 
active region. However, further study is needed to 
precisely quantify scaling properties associated to spe- 
cific active region features (e.g. emerging magnetic 
flux along the main polarity inversion line, sunspots 
build-up, delta-configuration...) and how the WTMM 
method can be sensitive to these elements. More pre- 
cisely, when analyzing higher resolution images than 
MDI, we expect that this segmentation tool will be all 
the more necessary as quiet-Sun and active region fea- 
tures are more entangled. 

The main outcome of the present study is the 
demonstration that the proposed multi-scale segmenta- 
tion procedure provides an objective way of studying 
the complexity in active regions separately from the 
surrounding quiet Sun. As such our results are signif- 
icantly more stable and robust when co mpared to pre 



vious fractal and mu l tifract al analysis (IPontieri et al. 
20031: iMcAteer etall l2005h . In a forthcoming pa- 



per, we will report on the application of this seg- 
mentation method to characterize the evolution of 
active regions keeping track of the multifractal pa- 
rameters f or possible correlatio ns with extreme so- 
lar events (iGallagher et al.l 120071) . Other applications 
of this method are in progress as the analysis of the 
intrinsic multifractal properties of entangled hot and 
cold int erstellar atomic gas from 3D numerical simu- 
lations dKestener & Auditll2009l) . 
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Fig. 12. — Multifractal analysis of a set of 5 active 
region magnetogram images (505x505) correspond- 
ing to the two sub- skeletons identified in Figures [TO] 
and[TT] The symbols (•) correspond to the segmented 
quiet Sun and (o) to the segmented active region, (a) 
h(q, a) vs log 2 a for different values of q\ the solid 
lines are linear regression fits over the range of scales 
a e [2°,2 30 ] cr w . (b) D(q,a) vs log 2 a. (c) r(q) vs 
q. (d) D(h) vs h. Error bars correspond to standard 
deviations in the linear regression procedure. 
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